Assessing crab bioturbation
First report: general overview and results
Importance
Bioturbation is the disturbance of sediments by biogenic reworking, and biodeposition is the deposition of biogenic material in the sediment.
These two processes are important in the soil and benthos dynamics because they increase ground surface area and particle mixing, which in turns enable and enhance ground oxygenation, gas interchange and nutrient recycling.
Crab bioturbation has been found to influence the vegetation structure, plant productivity and the bacterial community composition in the sediment.
Background
Traditionally, bioturbation in burrowing crabs is studied by analyzing burrow casts or the amount of loose sediment in the surface around a burrow. Both methods assume that bioturbation and biodeposition are a function of sediment removed and deposited elsewhere, and sediment compactation.
The burrow casts method uses the total volume of the burrow chamber and its surface area as a surrogate of bioturbation. Estimates of bioturbation rates (i.e. mass/volume per unit time) using this approach are difficult and prone to speculation as the absolute volume of the burrow represents the sediment moved and compacted over an undertermined number of days. Given its invasive and costly nature, this method is only applied to few burrows. Thus, this methods calculate bioturbation rates as the product of the total number of burrows and the average volume per burrow. Variable burrow diameter and architecture can introduce additional uncertainties in this method. The study of burrow casts is worth persuing, and it can inform and complement other methods, but by itself is likely to produce bioturbation and biodeposition estimates with large uncertainty intervals.
The amount of loose sediment is a direct proxy of the amount of sediment reworked by crabs. By using this method it is possible to obtain precise values of the amount of mass moved by crabs. This method requires frequent sampling, so loose sediment created by crab activities, can be identified and collected before is eroded or compacted. Collecting loose sediment in the surface around a burrow openning requires little equipment and intervation, so it can be done over relatively large areas. However, collecting the loose sediment deposited inside the burrow is more complicated (as far as I know it has not been done). Sampling the loose sediment inside the burrow is likely to affect the structurtal integrity of the chamber, and produce additional loose sediment which can increase sampling error. For this reason, some researchers (e.g. Wang et al 2010) have used artificial burrows (made of PVC) to measure the amount of loose sediment deposited inside burrows by crabs.
Collecting loose sediment is a better solution for assesing crab bioturbation. Sediment collected can also be used to answer additional questions about the consequences of bioturbation in the sediment properties. While this method require little equipment during sampling, it can be very time consuming when applied over large areas or areas with high burrow density.
The general assumption when studying crab bioturbation is that making a burrow requires excavating, compacting and/or depositing sediment in the burrow or elsewhere. Thus, the number of burrows and their size are believed to be related to bioturbation and biodeposition rates. In addition, crab feeding activity will also contribute to bioturbation and biodeposition.
In this paper, the method described in Herrera et al (2020)) is used to estimate bioturbation and evaluate potential processes underlying the observed bioturbation patterns.
Aims
-
Assess the bioturbation and biodeposition activity of burrowing crabs.
- Q: How much sediment matter (estimated as volume) do crabs move per unit time?
-
Study the potential process explaining bioturbation patterns.
-
Q: Do the number of burrows explain the observed change in volume?
-
Q: Do the distribution of burrows in the quadrat is associated with areas of major change in volume?
-
Q: Do the burrows distribution correlates with the local geometry of the sediment surface? This is, burrows position explain the observed changes in verticality, surface variation and other geometric features.
-
Methods
Here the computer vision and photogrammetry method described by (Herrera et al 2020) is used to assess the amount of change in the sediemnt due to burrowing and feeding activities by crabs.
Our method is a fundamental improvement to the estimation of crab bioturbation as it relates sediment turnover (i.e. moved matter) to observed volume changes in the sediment surface per unit time.
This study was performed in a mudflat next to the Ross River in Townsville, Queensland, Australia. This area was selected considering its easy access and the proven presence of fiddler crab from species Tubuca polita and Tubuca coarctata. The tidal regime in this area is well known and inundation can be predicted with precision.
Five permanent quadrats were installed in sediment areas exclusively or predominantly dominated by Tubuca polita individuals. Tubuca polita individuals are primarily active in the sediment surface during/following spring tides. To avoid confounding bioturbation with other natural process that can modify or moved matter, such as tidal erosion, the permanent quadrats were installed out of the area of tidal influence, this means, quadrats were not flooded.
All quadrats covered a 1 $m^{2}$ (1 x 1 m). Vertices from each quadrat were marked with Ground Control Points (GCPs). During the sampling time only two GCPs were lost. Replacements were installed inmediately after a lost was identified using a plastic reference quadrat.
Quadrats were followed during 14 days, and sampling occurred simultaneosly for all quadrats, same day within minutes of difference. First sampling event occurs at day zero (DO), with succesive samplings after one day (D1), three days (D3), seven days (D7), ten days (D10), and fourteen days (D14). Each sample event consisted in recording a video perpendicular to the sediment following the method described in Herrera et al 2020.
For each quadrat and sample event three products were created using Structure from Motion following Herrera et al (2020) methodos: point cloud data, textured and coloured mesh model and coloured orthomosaics.
Point Cloud Data
Point Cloud Data refers to coordinate data (x, y, z) created from Photogrammetry and other techniques. Point Cloud Data is the essential data type from which many 3D modelling and analysis tools based upon. This type of data represents a point in space based on its coordinates x, y, and z in a 3D space.
Rasters
As any other data set in higher dimensions (> 2D) statistical techniques can be used to reduce dimensionality. In many computer vision applications PCA is used to analyse point cloud data. For instance, to find features such as the edges of a geometry, or define the slope in an particular area.
Below there is a small representation of a point cloud data set. We could calculate a PCA for the entire dataset, and it will result in a new representation of the data in a 2D space. It is very likely that in such PCA the two principal components will be associated to axes "x" and "y", as much of the variation occurs in these two axes (i.e. there is not much changes in "z" or elevation). Thus, eigenvalues for "x" and "y" will be larger than "z" eigenvalues.
Alternatively we could recursively calculate PCAs in small windows including just a small subset of the data. This will produce many PCAs for small sections of the data. Areas where most variation occurs in the plane "xy" will show larger eigenvalues in "x" and "y". However, if the window is small enough, we could detect areas where change in elevation ("z") are more important than changes in the plane "xy". For instance, areas where the surface is wrinkled or corrugated, instead of flat. Therefore, in these areas the "z" eigenvalue will be larger.
Based on this principle, we can reduce our point data dimensionality to 2D (i.e. an image or raster). For example, an image could show for a given window or cell size the variation in the magnitude of the PCA1 as a colour gradient. In fact, and based on linear combinations of the eigenvalues, computer scientists have came with multiple metrics that describe geometric features of point cloud data. Some of these are: PCA1, PCA2, Planarity, Spherecity, Surface Variation, Verticality and 3rd Eigenvalue.
These features can be relevant when studying the effect of biological processes on the 3D space.
How do these features look like?
The geometric features described above are presented in the compound figure below for one video sample.
The first plot (top-left) shows the elevation (or height in cm) for this quadrat. Elevation represents the value of the "z" coordinate. The second plot (top-middle) shows roughness, which is a measure of the distance (cm) between a specific point and the plane made by the neighbors points inside a specific sphere window (radius = 1 cm). Essentially high roughness values indicate uneven or irrfular surfaces, i.e. points which are distance from the plane created by its neighbors points. All other plots represent geometric features that were calculated after calculating PCAs for a sphere cell/window of 1 cm radius. Features are arithmetic combinations from eigenvalues and eigenvectors (as per Hackel et al 2016).
Some geometric features are not very informative in this particular tasks. However, some can provide interesting insgights into the local shape/geometry in the surface. For instance, verticality (top-right) can be interpreted as changes in the sediment slope. Planarity, Surface Variation and Sphericity are also good descriptors of the sediment roughness.
#collapse-hide
options(repr.plot.width=10, repr.plot.height=10)
# Visualize burrows in quadrat and change in Digital Surface Models
## Load clean data
features_list <- readRDS(here::here('data/clean', 'features_list.rds'))
A_D0 <- features_list[[1]]
features <- names(A_D0)[c(4:length(names(A_D0)))]
plots <- list()
for (feat in features){
# Geometry feature plot
plots[[feat]] <- ggplot(A_D0, aes_string("x", "y", fill=feat)) +
geom_tile() +
scale_fill_viridis_c(name = paste0(feat)) +
coord_equal(xlim = c(0, 100), ylim = c(0, 100)) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
labs(x = feat) +
theme(# plot.margin = margin(5.5, 40, 5.5, 5.5),
legend.position = "right",
# legend.key.size = unit(0.8, "cm"),
legend.title = element_blank(),
axis.title.y = element_blank(),
# legend.justification = c("right", "top"),
panel.background = element_blank(),
axis.line = element_line(color = "black", size = .5),
panel.grid.major = element_line(color = "gray", size = 0.25),
legend.key = element_rect(fill = "white", color = NA))
}
plot_grid(plotlist = plots)
Textured and coloured mesh model
A second type of product that can be obtained from this workflow are textured and coloured mesh. In the current workflow used in this paper, these are only used for educational purposes (e.g. visualization for presentation). As most of the analysis are done using either point cloud data or the images created from point cloud data. The animation below shows one example of a textured and coloured mesh model.
Coloured orthomosaics
Coloured orthomosaics, a 2D image or raster, are created by combaning multiple information (i.e. pixels) from the video scan. This images were used to count the number of burrows within each quadrat at each time.
Analysis and results
Aim 1: Assess the bioturbation and biodeposition activity of burrowing crabs
- Question: How much sediment matter (estimated as volume) do crabs move per unit time?
To answer this question point cloud data was scaled and ortorectified. For each possible temporal combination within each quadrat (i.e. D0 vs D1, D0 vs D14, D1 vs D14) the amount of volumen change was calculated. This was done by aligned and superimposed pair of clouds using GCPs as coordinate references. Finally, the distance between the cloud pair was computed and the volume was estimated.
Three comparisons per quadrat were possible: D0 was compared to D1 and D14, and D1 was compared to D14. These comparisons resulted in volume change estimates for three time periods: one, thirteen and fourteen days.
The figure below shows the amount of volume added or removed for each quadrat (A to E) and time elapsed (1, 13, and 14 days). For most quadrats and time intervals the amount of added volume was larger than the amount of volume removed. The magnitude of this change was quite variable ranging from 143 to 8954 cubic centrimetres.
#collapse-hide
options(repr.plot.width=10, repr.plot.height=10)
# Visualize burrows in quadrat and change in Digital Surface Models
## Load clean data
biotur <- read.csv(here::here('data/raw', 'bioturbation.csv'), header = TRUE, sep = ",",
colClasses = c("factor", "factor", "numeric", "integer", "numeric", "numeric",
"numeric", "numeric", "numeric", "factor", "integer",
"numeric", "numeric", "numeric"))
biotur$quadrat <- as.factor(as.character(stringr::str_sub(biotur$Cloud.1, start = 5, end = 5)))
biotur$time1 <- as.factor(as.character(stringr::str_sub(biotur$Cloud.1, -7, -1)))
biotur$time2 <- as.factor(as.character(stringr::str_sub(biotur$Cloud.2, -7, -1)))
biotur$time1 <- sapply(biotur$time1, function(x) gsub("23NOV19", "D0", x), USE.NAMES = TRUE)
biotur$time1 <- sapply(biotur$time1, function(x) gsub("24NOV19", "D1", x), USE.NAMES = TRUE)
biotur$time1 <- sapply(biotur$time1, function(x) gsub("07DEC19", "D14", x), USE.NAMES = TRUE)
biotur$time2 <- sapply(biotur$time2, function(x) gsub("23NOV19", "D0", x), USE.NAMES = TRUE)
biotur$time2 <- sapply(biotur$time2, function(x) gsub("24NOV19", "D1", x), USE.NAMES = TRUE)
biotur$time2 <- sapply(biotur$time2, function(x) gsub("07DEC19", "D14", x), USE.NAMES = TRUE)
biotur$Bio_rate <- as.numeric(biotur$Bio_rate)
biotur$Time <- as.factor(as.integer(biotur$Time))
library(reshape2)
keep_names <- names(biotur)
keep_names <- keep_names[c(-5, -6)]
# keep_names
biotur1 <- melt(biotur, id.vars = keep_names)
biotur1$value <- ifelse(biotur1$variable == "Removed.volume", biotur1$value * -1, biotur1$value)
p1 <- ggplot(biotur1, aes(x = Time, y = value, fill=variable)) +
geom_bar(stat = "identity") +
facet_grid(~quadrat) +
scale_fill_manual("Volume", values = c("#848484", "#000000"), labels = c("Added", "Removed")) +
scale_y_continuous(breaks = seq(-3000, 9000, 1500)) +
labs(x = 'Time lapsed (days)', y = 'Change of volume in the sediment (cm^3)') +
theme(
axis.text = element_text(size = 15),
axis.title = element_text(size = 18),
# legend.position = "none",
panel.border = element_blank(),
panel.background = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
axis.line = element_line(color = "black", size = .5),
# legend.key = element_rect(fill = "white", color = NA),
axis.text.x = element_text(hjust = 1),
strip.placement = "inside",
strip.text = element_text(size = 18),
strip.background = element_blank(),
legend.text = element_text(size = 18),
legend.title = element_text(size = 20))
p1
The bioturbation rate was calculated as the absolute value resulting from the arithmetic sum of the added and removed volume per square meter per day. The bioturbation rate shows that this estimate is quite sensitive to the time scale used. This means that while the net amount of volume change after 13 or 14 days is quite high (figure above), when values are standardized by unit time, the resulting rate estimates ended being quite small.
This result is explained (or confounded) with the observed pattern of decreasing surface activity of fiddler crabs towards the end of the experiment. Tubuca polita fiddler crabs were rarely observed in the surface after after six days since the experiment initial samplig. This species is mostly active around the spring tides period.
#collapse-hide
options(repr.plot.width=8, repr.plot.height=5)
p2 <- ggplot(biotur1, aes(x = Time, y = abs(Bio_rate_1sqm))) +
geom_point(size = 5) +
facet_grid(~quadrat) +
scale_x_discrete("Time interval considered (days)") +
scale_y_continuous("Bioturbation rate (cm3 per m2d-1)", breaks = seq(500, 3500, 500)) +
theme(
axis.text = element_text(size = 15),
axis.title = element_text(size = 18),
# legend.position = "none",
panel.border = element_blank(),
panel.background = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
axis.line = element_line(color = "black", size = .5),
# legend.key = element_rect(fill = "white", color = NA),
axis.text.x = element_text(angle = 90, size = 14),
strip.placement = "inside",
strip.text = element_text(size = 18),
strip.background = element_blank(),
legend.text = element_text(size = 18),
legend.title = element_text(size = 20))
p2
The animation below shows the change in the sediment in one quadrat after one day.
Aim 2: Study the potential process explaining bioturbation patterns
- Question: Do the number of burrows explain the observed change in volume?
The total number of burrows decreased through time in all quadrats. This was expected after the observed reduction in the activity of crabs in the field through time. Quadrats B and C showed noticeble less burrows than A, C and D.
#collapse-hide
options(repr.plot.width=8, repr.plot.height=6)
burrows <- readRDS(here::here('data/clean', 'burrows_master.rds'))
burrows$quadrat <- as.factor(as.character(stringr::str_sub(burrows$ID, 15, 15)))
burrows$day <- as.factor(as.character(stringr::str_sub(burrows$ID, -3, -1)))
burrows$day <- gsub("_", "", burrows$day)
# Plot total number of burrows per quadrat and time
burrow_tally <-burrows %>%
dplyr::group_by(quadrat, day) %>%
dplyr::tally()
b_tally_plot<- ggplot(burrow_tally, aes(day, n)) +
geom_bar(position = "dodge", stat = "identity", fill = "black") +
facet_wrap(~ quadrat, scales = "free_x", nrow = 1) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 250), breaks = seq(0, 240, 20)) +
# coord_flip() +
labs(x = 'Time', y = 'Number of burrows') +
theme(axis.text = element_text(size = 10),
legend.position = "none",
panel.border = element_blank(),
panel.background = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
axis.line = element_line(color = "black", size = .5),
legend.key = element_rect(fill = "white", color = NA))
b_tally_plot
As the estimates of bioturbation rate are sensitive to the time interval used, the figures below show the relationship between the number of burrows and the bioturbation rate calculated using different time intervals. This relationship was explored using the non-transformed and log transformed bioturbation rate.
The first pair of figures show this relationship using the number of burrows at the beginning of the time interval as predictor of bioturbation rate. There is not a robust pattern. The bioturbation rate calculated after one day seems to be positive related to the number of burrows.
#collapse-hide
options(repr.plot.width=12, repr.plot.height=6)
## Load clean data
## Load clean data
burrows <- readRDS(here::here('data/clean', 'burrows_master.rds'))
library(dplyr)
burrows_tally <- burrows %>%
group_by(ID) %>%
tally() %>%
mutate(quadrat = as.factor(as.character(str_remove_all(ID, "burrowmap_C01_")))) %>%
mutate(quadrat = as.factor(as.character(stringr::str_sub(quadrat, 1, 1)))) %>%
mutate(quadrat_time = as.factor(as.character(stringr::str_sub(ID, -3, -1)))) %>%
mutate(quadrat_time = gsub("_", "", quadrat_time))
# burrows_tally$time <- rep(seq(1,3,1), 5)
# biotur1$nburrows1 <- burrows_tally$n[match(paste(biotur1$quadrat, biotur1$time1),
# paste(burrows_tally$quadrat, burrows_tally$quadrat_time)
# )]
biotur1$nburrows1 <- burrows_tally$n[match(interaction(biotur1$quadrat, biotur1$time1),
interaction(burrows_tally$quadrat, burrows_tally$quadrat_time))]
biotur1$nburrows2 <- burrows_tally$n[match(interaction(biotur1$quadrat, biotur1$time2),
interaction(burrows_tally$quadrat, burrows_tally$quadrat_time))]
# biotur1$nburrows1 <- ifelse(is.na(biotur1$nburrows1), biotur1$nburrows2, biotur1$nburrows1)
biotur1$burrow_rate <- round(abs((biotur1$nburrows2-biotur1$nburrows1) / biotur1$nburrows1), digits = 2)
biotur1$comparison <- paste0(biotur1$time1, "_vs_", biotur1$time2)
# Observe that I use the variable value "Removed.volume", "DX" to subset the data and obtain the right pair combinations,
# but I am actually plotting bio_rate_1sqm
p3 <- ggplot() +
geom_point(data = subset(biotur1, variable == "Removed.volume" & time2 == "D1"),
aes(x = nburrows1, y = log(abs(Bio_rate_1sqm)), colour = "col1"), size = 5) +
geom_point(data = subset(biotur1, variable == "Removed.volume" & time1 == "D1" & time2 == "D14"),
aes(x = nburrows1, y = log(abs(Bio_rate_1sqm)), colour = "col2"), size = 5) +
geom_point(data = subset(biotur1, variable == "Removed.volume" & time1 == "D0" & time2 == "D14"),
aes(x = nburrows1, y = log(abs(Bio_rate_1sqm)), colour = "col3"), size = 5) +
scale_colour_manual("Time elapsed", labels = c("1 day", "13 days", "14 days"),
values = c("col1" = "#bf812d", "col2" = "#35978f", "col3" = "#000000")) +
labs(x = "Number of burrows at the beginning time", y = "Bioturbation rate ( log(cm3 per m2d-1) )") +
theme(axis.text = element_text(size = 15),
axis.title = element_text(size = 18),
legend.position = "top",
panel.border = element_blank(),
panel.background = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
axis.line = element_line(color = "black", size = .5),
# legend.key = element_rect(fill = "white", color = NA),
axis.text.x = element_text(size = 14),
# strip.placement = "inside",
# strip.text = element_text(size = 18),
# strip.background = element_blank(),
legend.text = element_text(size = 18),
legend.title = element_text(size = 20))
# Observe that I use the variable value "Removed.volume", "DX" to subset the data and obtain the right pair combinations,
# but I am actually plotting bio_rate_1sqm
p4 <- ggplot() +
geom_point(data = subset(biotur1, variable == "Removed.volume" & time2 == "D1"),
aes(x = nburrows1, y = abs(Bio_rate_1sqm), colour = "col1"), size = 5) +
geom_point(data = subset(biotur1, variable == "Removed.volume" & time1 == "D1" & time2 == "D14"),
aes(x = nburrows1, y = abs(Bio_rate_1sqm), colour = "col2"), size = 5) +
geom_point(data = subset(biotur1, variable == "Removed.volume" & time1 == "D0" & time2 == "D14"),
aes(x = nburrows1, y = abs(Bio_rate_1sqm), colour = "col3"), size = 5) +
scale_colour_manual("Time elapsed", labels = c("1 day", "13 days", "14 days"),
values = c("col1" = "#bf812d", "col2" = "#35978f", "col3" = "#000000")) +
labs(x = "Number of burrows at the beginning time", y = "Bioturbation rate (cm3 per m2d-1)") +
theme(axis.text = element_text(size = 15),
axis.title = element_text(size = 18),
legend.position = "top",
panel.border = element_blank(),
panel.background = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
axis.line = element_line(color = "black", size = .5),
# legend.key = element_rect(fill = "white", color = NA),
axis.text.x = element_text(size = 14),
# strip.placement = "inside",
# strip.text = element_text(size = 18),
# strip.background = element_blank(),
legend.text = element_text(size = 18),
legend.title = element_text(size = 20))
plot_grid(p3, p4)
The next two figures use the number of burrows at the end of the interval as predictor of bioturbation rate. Again there is not a robust pattern. But again the bioturbation rate calculated after one day seems to be related to the number of burrows.
#collapse-hide
options(repr.plot.width=12, repr.plot.height=6)
# Observe that I use the variable value "Removed.volume", "DX" to subset the data and obtain the right pair combinations,
# but I am actually plotting bio_rate_1sqm
p5 <- ggplot() +
geom_point(data = subset(biotur1, variable == "Removed.volume" & time2 == "D1"),
aes(x = nburrows2, y = log(abs(Bio_rate_1sqm)), colour = "col1"), size = 5) +
geom_point(data = subset(biotur1, variable == "Removed.volume" & time1 == "D1" & time2 == "D14"),
aes(x = nburrows2, y = log(abs(Bio_rate_1sqm)), colour = "col2"), size = 5) +
geom_point(data = subset(biotur1, variable == "Removed.volume" & time1 == "D0" & time2 == "D14"),
aes(x = nburrows2, y = log(abs(Bio_rate_1sqm)), colour = "col3"), size = 5) +
scale_colour_manual("Time elapsed", labels = c("1 day", "13 days", "14 days"),
values = c("col1" = "#bf812d", "col2" = "#35978f", "col3" = "#000000")) +
labs(x = "Number of burrows at the end time", y = "Bioturbation rate ( log(cm3 per m2d-1) )") +
theme(axis.text = element_text(size = 15),
axis.title = element_text(size = 18),
legend.position = "top",
panel.border = element_blank(),
panel.background = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
axis.line = element_line(color = "black", size = .5),
# legend.key = element_rect(fill = "white", color = NA),
axis.text.x = element_text(size = 14),
# strip.placement = "inside",
# strip.text = element_text(size = 18),
# strip.background = element_blank(),
legend.text = element_text(size = 18),
legend.title = element_text(size = 20))
# Observe that I use the variable value "Removed.volume", "DX" to subset the data and obtain the right pair combinations,
# but I am actually plotting bio_rate_1sqm
p6 <- ggplot() +
geom_point(data = subset(biotur1, variable == "Removed.volume" & time2 == "D1"),
aes(x = nburrows2, y = abs(Bio_rate_1sqm), colour = "col1"), size = 5) +
geom_point(data = subset(biotur1, variable == "Removed.volume" & time1 == "D1" & time2 == "D14"),
aes(x = nburrows2, y = abs(Bio_rate_1sqm), colour = "col2"), size = 5) +
geom_point(data = subset(biotur1, variable == "Removed.volume" & time1 == "D0" & time2 == "D14"),
aes(x = nburrows2, y = abs(Bio_rate_1sqm), colour = "col3"), size = 5) +
scale_colour_manual("Time elapsed", labels = c("1 day", "13 days", "14 days"),
values = c("col1" = "#bf812d", "col2" = "#35978f", "col3" = "#000000")) +
labs(x = "Number of burrows at the end time", y = "Bioturbation rate (cm3 per m2d-1)") +
theme(axis.text = element_text(size = 15),
axis.title = element_text(size = 18),
legend.position = "top",
panel.border = element_blank(),
panel.background = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
axis.line = element_line(color = "black", size = .5),
# legend.key = element_rect(fill = "white", color = NA),
axis.text.x = element_text(size = 14),
# strip.placement = "inside",
# strip.text = element_text(size = 18),
# strip.background = element_blank(),
legend.text = element_text(size = 18),
legend.title = element_text(size = 20))
plot_grid(p5, p6)
I have also used the change in number of burrows as a predictor of bioturbation rate. The proportional change in number of burrows was calculated as one minus the ratio between the final number of burros and initial number of burrows. There is not an evident relationship.
#collapse-hide
options(repr.plot.width=16, repr.plot.height=8)
# Observe that I use the variable value "Removed.volume", "DX" to subset the data and obtain the right pair combinations,
# but I am actually plotting bio_rate_1sqm
p7 <- ggplot() +
geom_point(data = subset(biotur1, variable == "Removed.volume" & time2 == "D1"),
aes(x = burrow_rate, y = log(abs(Bio_rate_1sqm)), colour = "col1"), size = 5) +
geom_point(data = subset(biotur1, variable == "Removed.volume" & time1 == "D1" & time2 == "D14"),
aes(x = burrow_rate, y = log(abs(Bio_rate_1sqm)), colour = "col2"), size = 5) +
geom_point(data = subset(biotur1, variable == "Removed.volume" & time1 == "D0" & time2 == "D14"),
aes(x = burrow_rate, y = log(abs(Bio_rate_1sqm)), colour = "col3"), size = 5) +
scale_colour_manual("Time elapsed", labels = c("1 day", "13 days", "14 days"),
values = c("col1" = "#bf812d", "col2" = "#35978f", "col3" = "#000000")) +
labs(x = "Proportion of burrows left at the end of interval", y = "Bioturbation rate ( log(cm3 per m2d-1) )") +
theme(axis.text = element_text(size = 15),
axis.title = element_text(size = 18),
legend.position = "top",
panel.border = element_blank(),
panel.background = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
axis.line = element_line(color = "black", size = .5),
# legend.key = element_rect(fill = "white", color = NA),
axis.text.x = element_text(size = 14),
# strip.placement = "inside",
# strip.text = element_text(size = 18),
# strip.background = element_blank(),
legend.text = element_text(size = 18),
legend.title = element_text(size = 20))
# Observe that I use the variable value "Removed.volume", "DX" to subset the data and obtain the right pair combinations,
# but I am actually plotting bio_rate_1sqm
p8 <- ggplot() +
geom_point(data = subset(biotur1, variable == "Removed.volume" & time2 == "D1"),
aes(x = burrow_rate, y = abs(Bio_rate_1sqm), colour = "col1"), size = 5) +
geom_point(data = subset(biotur1, variable == "Removed.volume" & time1 == "D1" & time2 == "D14"),
aes(x = burrow_rate, y = abs(Bio_rate_1sqm), colour = "col2"), size = 5) +
geom_point(data = subset(biotur1, variable == "Removed.volume" & time1 == "D0" & time2 == "D14"),
aes(x = burrow_rate, y = abs(Bio_rate_1sqm), colour = "col3"), size = 5) +
scale_colour_manual("Time elapsed", labels = c("1 day", "13 days", "14 days"),
values = c("col1" = "#bf812d", "col2" = "#35978f", "col3" = "#000000")) +
labs(x = "Proportion of burrows left at the end of interval", y = "Bioturbation rate (cm3 per m2d-1)") +
theme(axis.text = element_text(size = 15),
axis.title = element_text(size = 18),
legend.position = "top",
panel.border = element_blank(),
panel.background = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
axis.line = element_line(color = "black", size = .5),
# legend.key = element_rect(fill = "white", color = NA),
axis.text.x = element_text(size = 14),
# strip.placement = "inside",
# strip.text = element_text(size = 18),
# strip.background = element_blank(),
legend.text = element_text(size = 18),
legend.title = element_text(size = 20))
plot_grid(p7, p8)
- Question: Do the distribution of burrows in the quadrat is associated with areas of major change in volume?
The first step to analyse this question is to evaluate if the burrow distribution follows what can be described as a complete spatial randomness process. This is, are burrows randomly distributed in the quadrat, or alternitavely these are ordely scattered, or clustered.
The number of burrows inside a quadrat can be thought as a random Poisson variable. Therefore, we can use a Poisson Point Process as a the null model (i.e. complete spatial randomness). This model assumes that a burrow can be located in any position inside the quadrat and that the position of a burrow is not affected by the existence and position of other burrows.
The two graphs below show the same number of burrows, but one was simulated using a Poisson point process (complete spatial randomness, left), and the other one was created using the observed burrow distribution for one sample (right). It is important to notice that often, as in the case of the two images below, is difficult to visually differentiate between a Poisson point process, i.e. complete spatial randomness, and any other process creating a clustered burrow distribution.
#collapse-hide
options(repr.plot.width=8, repr.plot.height=4)
# Visualize burrows distribution in quadrats
## Load clean data
burrows <- readRDS(here::here('data/clean', 'burrows_master.rds'))
## Create factors for quadrat and time
burrows$quadrat <- as.factor(as.character(stringr::str_sub(burrows$ID, 15, 15)))
burrows$day <- as.factor(as.character(stringr::str_sub(burrows$ID, -3, -1)))
burrows$day <- gsub("_", "", burrows$day)
# Generate a random point pattern
set.seed(1024)
poi_point <- runifpoint(145)
poi_point <- data.frame(x = poi_point$x*100, y = poi_point$y*100)
pa <- ggplot(poi_point, aes(x, y)) +
geom_point(size = 2) +
coord_equal(xlim = c(0, 100), ylim = c(0, 100)) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
labs(x = 'Distance (cm)', y = 'Distance (cm)') +
theme(# plot.margin = margin(5.5, 40, 5.5, 5.5),
legend.position = "right",
panel.border = element_rect(fill = NA),
panel.background = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
axis.line = element_line(color = "black", size = .5),
legend.key = element_rect(fill = "white", color = NA))
pb <- ggplot(subset(burrows, quadrat == "A" & day == "D0"), aes(x, y)) +
geom_point(size = 2) +
coord_equal(xlim = c(0, 100), ylim = c(0, 100)) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
labs(x = 'Distance (cm)', y = 'Distance (cm)') +
theme(# plot.margin = margin(5.5, 40, 5.5, 5.5),
legend.position = "right",
panel.border = element_rect(fill = NA),
panel.background = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
axis.line = element_line(color = "black", size = .5),
legend.key = element_rect(fill = "white", color = NA))
plot_grid(pa, pb)
The burrow distribution for all quadrats and times (next figure) shows a reduction in the number of burrows through time, as already observed in a previous figure. Some quadrats seem to contain cluster of burrows.
#collapse-hide
# options(repr.plot.width=16, repr.plot.height=8)
burrows_map <- ggplot(burrows, aes(x, y)) +
geom_point(size = 0.8) +
facet_grid(day ~ quadrat) +
coord_equal(xlim = c(0, 100), ylim = c(0, 100)) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
labs(x = 'Distance (cm)', y = 'Distance (cm)') +
theme(# plot.margin = margin(5.5, 40, 5.5, 5.5),
legend.position = "right",
panel.border = element_rect(fill = NA),
panel.background = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
axis.line = element_line(color = "black", size = .5),
legend.key = element_rect(fill = "white", color = NA))
burrows_map <- burrows_map + theme(panel.spacing = unit(10, "points"))
burrows_map
If burrows follow a stationary Poisson process, this is, burrows are equally likely to be found in any space in the quadrat, then we would expect that the probability of fiding two burrows at any given distance is the same. This probability can be modelled for the null hypothesis (i.e. stationary Poisson process) and it can be calculated for the observed data. One method to estimate if the observed data deviate from a stationary Poisson process is the Pair Correlation Function (Stoyan & Stoyan 1994).
The Pair Correlation Function (PCF) was calculated for all samples and considering distances among burrows from 0 to 25 centimeters. In the figure below, the x axis represent all distances among burrows that were evaluated. The pair correlation fucntion, y axis, represents the normalized probability of finding two burrows at locations a and b seperated by distance "r". The horizontal dotted red line crossing the y axis at 1 represents the expected PCF value for a complete random Poisson process. The red envelope around this line is the calculated significance band. The continuous black line represents the observed probability of finding two burrows at locations a and b.
In general g(r) values smaller than one indicate that the distance among burrows for any given distance "r" are less common than what would be expected by a complete random process. PCF ( g(r) ) values greater than one indicate the opposite, which suggest that burrows are clustered.
For most samples, the PCF suggest that burrows are more clustered than what would be expected by a random process in distances of 2-3 centimeters.
#collapse-hide
options(repr.plot.width=10, repr.plot.height=16)
## Load clean data
pcf.env_list <- readRDS(here::here('data/clean', '08_pcf.env_list.rds'))
pcf_names <- names(pcf.env_list)
pcf_plot_list <- c()
for (name_item in pcf_names){
b <- pcf.env_list[[name_item]]
b_plot <- ggplot(b) +
geom_line(aes(x = r, y = obs, colour = "black")) +
geom_line(aes(x = r, y = theo, colour = "red"), linetype = "dotdash") +
geom_ribbon(aes(x = r, ymin = lo, ymax = hi), fill = "red", alpha = 0.25) +
labs(x = "Radii (cm)", y = expression(paste("Pair correlation function ", italic("g(r)")))) +
scale_x_continuous(expand = c(0, 0), breaks = seq(0, 30, 5)) +
scale_y_continuous(breaks = c(0, 1, seq(5, max(b$hi, na.rm = TRUE), 5))) +
scale_color_manual("", labels = c(expression(paste(hat(g)["obs"], "(r)")),
expression(paste(hat(g)["theoretical"], "(r)"))), values = c("black", "red")) +
ggtitle(name_item) +
theme(
legend.text = element_text(size = 10),
legend.text.align = 0,
legend.position = c(0.8, 0.8),
panel.background = element_blank(),
axis.line = element_line(color = "black", size = .5),
panel.grid.major = element_blank(),
legend.key = element_rect(fill = "white", color = NA),
plot.title = element_text(size = 15, hjust = 0.2, vjust = -10))
pcf_plot_list[[name_item]] <- b_plot
rm(list = c("name_item", "b", "b_plot"))
gc()
}
plot_grid(plotlist = pcf_plot_list, ncol = 3)
The next step is to visualize how the potential point process underlying the burrow distribution would look like.
and then try to fit a point process model to the observed point pattern using covariates.
Importantly the assumption of homogeneity should be checked.
- Question: Do the burrows distribution correlates with the local geometry of the sediment surface? This is, burrows position explain the observed changes in verticality, surface variation and other geometric features.
References
Diggle, Peter J. Statistical analysis of spatial and spatio-temporal point patterns. CRC press, 2013.
Baddeley, Adrian, et al. Spatial Point Patterns : Methodology and Applications with R, CRC Press LLC, 2015. ProQuest Ebook Central, https://ebookcentral.proquest.com/lib/jcu/detail.action?docID=4742712.